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various integration schemes on the element performance are presented. Convergence 


studies for laminated composites for different fiber orientations are also discussed. 
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A bi-quadratic isoparametric plate/shell bending finite element is developed 
to study the behavior of isotropic and laminated composite plates. The element 
is based on Mindlin-Reissner's theory and the principle of virtual displacements. 
The element is implemented in a computer program. Results are presented and 
compared with analytical solutions to validate this element. Good agreement is 
observed for thin plates, while discrepancies are noted for thick plates. Effects of 
various integration schemes on the element performance are presented. Convergence 


studies for laminated composites for different fiber orientations are also discussed. 





I. INTRODUCTION 


A. THE SCOPE OF THE THESIS 
Finite element analysis provides a general tool to solve problems in struc- 
tural mechanics. The methodology is applicable for static and dynamic response of 
structures and in predicting the elastic stability limits. 
The focus of the present study is to develop tools to analyze laminated com- 
posite plates and validate the model by comparing with known solutions. 
More specifically, the objectives of the present study are: 
a) to review some of the pertinent literature in the area of laminated 
composite plates. 
b) develop a finite element for the analysis of composite plates. 
c) develop consistent mass and load matrices. 
d) study the effect of thickness to characteristic length of the plate. 
e) study the effects of integration schemes. 
The outline of the remainder of this thesis is as follows: 
The basic formulation of the stiffness, mass and load matrix for the bending 
of flat plates using laminated composite materials are described in Chapter II. 
Chapter III addresses certain aspects of computational implementation. 
Chapter IV describes some test cases, example calculations and comparison 
with classical plate theory. 
Finally, Chapter V reflects experience gained and some suggestions for future 


research. 


B. LITERATURE SURVEY 

The finite element method, [Ref. 1], may be described as a general discretiza- 
tion procedure of continuum problems, posed by mathematically defined statements 
with applications to several engineering analysis problems. 

A brief literature review pertaining to the analysis of plates/shells using finite 
element approximation is presented in the following paragraphs. A significant con- 
tribution to include shear is given by Mindlin (Ref. 2], while Hughes etal, [Ref. 3] 
adapt this theory to develop finite elements for the analysis of isotropic plates [Refs. 
3, 4]. Laminated plate theory based on the classical Kirchoff hypothesis has been 
established by Reissner and Starsky [Ref. 5] and Whitney and Leissa [Ref. 6]. The 
effect of reduced integration in isoparametric elements was presented by Zienkiewicz 
etal [Ref. 7] and Hughes etal [Ref. 3]. 

The finite element method of analysis for the plate bending problem including 
shear deformation has been presented by Pryor and Barker [Ref. 8]. Mawenya 
describes formulations for multi-layer plates [Refs. 9, 10]. The higher order shear 
deformation theory of laminated composite plates was developed by Krishna Murty 
(Ref. 11], and Lo et.al. [Refs. 12, 13] present a higher order, three-dimensional 
theory. 

Burt [Ref. 14] presented a higher order theory and compared with Pagano's 
elasticity-theory solution for the case of cylindrical bending and a symmetric cross- 
ply laminate consisting of three equal-thickness layers. Bending of simply supported 
thick rectangular plates was presented by Srinivas and Rao [Ref. 15]. Exact elas- 
ticity solutions for some particular plate bending problems have been obtained by 


Pagano [Refs. 16, 17, 18, 19] and Srinivas and Rao [Ref. 20]. 


Application of classical shell theory, including transverse shear deformation is 
presented by Vinson and Chou [Ref. 21]. Naschie [Ref. 22] studied large deflec- 
tion behavior of orthotropic composite materials. Srinivas [Ref. 23] developed a 
refined approximate theory for the static and dynamic analysis of finite, laminated, 
composite, circular cylindrical shells with general boundary conditions. 

Plate theories, which include shear deformation has been developed by Whit- 
ney [Ref. 24] and Mau [Ref. 25]. The first such theory for laminated isotropic 
plates is due to Yang, Norris and Stravinsky [Ref. 26]. Reddy [Ref. 27] developed 
a higher order shear deformation theory of laminated composite plates. Reddy and 
Sandidge [Ref. 28] presented mixed finite element models of the classical and shear 
deformation theory. The effect of transverse shear deformation on bending of elastic 
symmetric laminated composite plate undergoing large deformation is presented by 
Gorji [Ref. 29]. 

Based on anisotropic elasticity, Hearmon [Ref. 30] and Lekhnitskii [Ref. 31] 
present general theories for laminates. The covariant form for the transformed lam- 
ina stiffnesses has been given by Tsai and Pagano [Ref. 33]. Gibert and Schneider 
[Ref. 34] directed their study in this direction. Noor and Mathers [Refs. 35, 36] pre- 
sented the effects of shear deformation and anisotropy on the response of laminated 
anisotropic plates. 

Nelson and Lorch [Ref. 40] compare the accuracy of various plate models to 


predict the behavior of laminated orthotropic plates. 


II. THEORETICAL FORMULATION 


A. INTRODUCTION 


In this chapter, the derivation of plate finite elements based on Mindlin's 
theory is described. The principle of virtual displacements is invoked to obtain 


equilibrium relations. 


B. THE PRINCIPLE OF VIRTUAL WORK 


In this section, we prove that total internal virtual work is equal to total ex- 
ternal virtual work and equivalence of this principle to the minimum total potential 
energy principle. 

In general, the total potential energy of a structural system is equal to the 


sum of strain energy and potential energy. 


Il, -U4V (2.1) 


where, 
Il; : Total Potential Energy of Structural System 
U  : Total Strain Energy 
V Total Potential of External Loads 
The total minimum potential energy requires that first variation of total po- 


tential energy be zero, or 


ST, =0 (2.2) 


or, 


U +0V=0 (2.3) 


in other words, 


6U = —6V (2.4) 


This may also be written in a different form, recognizing U as being the work 
done by internal forces and that work done by the external forces being equal to the 


negative of the total potential energy of the external loads. That is, 


Win: = Wert (2.5) 


which is a statement of the principle of virtual work, stating that, if the body is in 
equilibrium, the total virtual work done by the internal forces is equal to the total 
virtual work done by external forces for arbitrary, kinematically admissible virtual 
displacements. 

It may be noted that the form in Equation (2.4) is restricted to conservative 
loadings while the form in Equation (2.5) is applicable for any loading form. 


The total internal virtual work may be written as 


Win = | to)" {5e}d(A) (2.6) 


where, 


fo) : Vector of Stresses 
{de} : Vector of Virtual Strains 


By using generalized Hooke’s law for material constitutive relations, stresses 


may be expressed as 


10) = [Qlie) (2.7) 
where the matrix [Q] contains the material stiffness coefficients. 


If the thickness t is constant, then total internal virtual work takes the form, 


me I t{e}T [Q] {5e}d(vol) (2.8) 

In order to derive the element matrices, the principle of virtual work, which is 

equally applicable to the element as well as to the total structure, is applied to the 

element. The virtual work is additive and results in the virtual work of the entire 
structure under consideration. 


The linear strain-displacement relations are given by 


{e} < [Bltu) (2.9) 


while the virtual strains are given by 


(če) — [B] (ču) (2.10) 


The operator matrix, [B], is dependent on the shape functions and their deriva- 
tives, and {u} and {ôu} are vectors of displacements and virtual displacements, 
respectively of the element. 


On substituting Equation (2.9), (2.10) in Equation (2.8), we obtain, 


SWine = | tóu} [BIQBH 6u} )dA (2.11) 

or, rearranging 
5Wine = {8u}"(t | [BF[QJIBJdA) (2.12) 
Win: = (6u)"(K] (2.13) 


where [K] = t [,[B"][Q][B]dA is the element stiffness matrix. 


In order to derive mass and load matrices, we start from the expression for 


external virtual work, 


sa f T46u)ds — / z(6u)d(vol) + {6u} {F} (2.14) 


where, T is the external force per unit length along the boundary of the element, 


X is the body force per unit volume (inertia force) 
{F} is the point loads applied at nodal points. 


On substituting for displacements in terms of nodal displacements using shape 


functions, we obtain, 


ŚW: = {6u}"( | [N]Tds) - 48u)"( | tpiNJINI dA)lu] 4 (8u) 4F) (245) 


By invoking the principle of virtual displacements and noting that {du} is 


arbitrary, we get the equilibrium equations in the following form; 


[K]{u} + [M]{u} = [F] + [F]" (2.16) 
where, 
[M] = T tp[N][N]TdA (2.17) 
F 
[F] =<: (2.18) 
Fs 
[ FJ? = / [N]T ds (2.19) 


It may be noted that the matrix [M] is the mass matrix, [F] is the vector of 


point loads and [F]" is the vector of consistent loads. 


C. LAMINATE THEORY 
1. Introduction 
In the expression for [K| matrix, there are three unknown matrices. In 
this section, we discuss about calculation of [Q] matrix. The matrix [Q] , relat- 
ing stresses and strains, consists of material stiffness coefficients. It reflects the 
properties of both fibers and matrix. 
The behavior of laminated plates is characterized by possible coupling 
between membrane action and bending action. 
Following discussions review some aspects of the laminate theory. 
2. Strain-displacement Relations 


As shown in Figure 2.1, a general strain field converts configuration 012 


10 0 | 2 . 
Linearized normal strains may be written as: 
Qu 
r = = 2.20 
e (2.20) 
and 
Qv 
CE — 2.21 
y Oy ( ) 


The shear strain is defined as the amount of change in a right angle. For 


small angles, 


dry = Pit Pr (2.22) 
Qu Ov 
= oF sp BT (223) 


Similar expressions may be derived for other strain components. 


The foregoing strain-displacement relations can be summarized in a matrix- 


operator form. 


— 0 0 
Ey 0 B 0 u 
z i2|l00 £|iv (2.24) 
ry 3 8 Q 1D 
y Oy Ox 
“| lo #2 
Ver a n 의 
02 Ox 


3. Stress-Strain Relations 


The generalized Hooke's Law takes the form, 


{o} = [Qlte) (2.25) 


where, the stress vector is given by 


CO aTa] (2.26) 


and the strain vector is given by 


E A a r] (221) 


The material stiffness matrix [Q] , contains nine independent coefficients. 
An orthotropic material has only five independent elastic coefficients — E, E2, v12, V21, and G12. 


More explicitly, Eguation (2.25) may be written as 


Or Qu Qu Gi 


0 0 0 Er 
Oy Qn Qo Q»5 0 0 0 €y 
O; = | Qa Qo Qs 0 0 0 Ez (2.28) 
Try 0 0 0 Q44 0 0 Yry 
Tyz 0 0 0 0 Qss 0 Yyz 
Ter 0 0 0 0 0 Ose Yzr 





Figure 2.1: Displacement and Distortion of Differential Lengths dx and 
dy. 


For planar orthrotropic material, the stress-strain relation reduces to: 


Or Qu (212 0 0 0 Ex 
Oy m Qa (Q22 0 0 0 Ey 2.29 
Ty 0 0 044 0 0 Dry SZ 
Tus U U U Qss 0 Yyz 
Ter J 0 0 0 0 0266 ] Ver 
where, 
Ei 
œ DE 
Qn = > (2.30) 
E, 
— AA 2.31 
Qn = pt (2.31) 
Qu = Gir (232) 
Qss = Gz (2:33) 
Qs = Gis (2.34) 


It may be noted that shears 7,, and 7,, are obtained to account for 

transverse shears. 
4. Lamina of Composite Materials 

Composite structures are built of individual lamina, which are stacked 
into several number of layers to form a laminate. Each lamina consists of, typically, 
uniaxial fibers embedded in a matrix, such as a resin. In Figure 2.2, the principal 
material axes are labelled 1 and 2, that is, 1-direction is parallel to the fibers direction 
and the 2-direction is normal to them. 

It may be noted that in each lamina there exists a state of plane stress. 

The state of stress is also shown in the Figure 2.2, in both 1-2 and x-y 
coordinate system. The computation stresses in different coordinate system follows 


the usual transformation rules [Ref. 21]. We follow the notation that 0,02 are 


11 


^ O, 





Figure 2.2: Lamina Coordinate System (2-Dimension) 


normal stresses, 03,04 and os are the shear stress in 1-2 system. €1, €2, €3,€4 and €5 


are the corresponding strains in 1-2 system. 


O) Or 
02 |=[Tlo, | vv (2.35) 
04 Try 


where the transformation matrix is given by 


m? n?  2mn 


[T],, =| n m? —2mn (2.36) 
-mn mn (m-n?) 


The direction cosines of the unit normal are determined from 


m = cos 0 (2:99) 


n = sinf (2.38) 


The subscript CL refers to two-dimensional case, that is the x-y plane 


only. Similarly, the strains are related in the two coordinate systems by 


61 Er 
a NE e, (2.39) 
ć12 Jzy 


By including transverse shears, we modify these transformations as fol- 


lows: 


O) Or 
02 Oy 
os | = [T] | Tay (2.40) 
04 Tyz 
05 Tzr 
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€2 

€3 = [T] Yry (2.4 l ) 
€4 Yyz 

€5 Jzz 


where the transformation is given by, 


m n*  2mn 0 0 
n? m? -2mn 0 0 
[I] 2| —=mn mn m*-n?*? 0 0 (2.42) 
0 0 0 m —n 
0 0 0 n m 


The stresses and strains in x-y system may simply be obtained from 1-2 


system by inverting the matrix [T]. 


Or 01 
O, 02 
Ty |=[T] < 03 (2.43) 
Tyz 04 
722 05 
and, 
Cr 61 
Ey €2 
Yay | = [IT] | es ps 
Yyz €4 
Yer . €5 
with 
m? mn? —2mn .0 0 
n? m? mn 0 0 
[27] < l mn —=mn m—-n 0 0 (2.45) 
0 0 0 m n 
0 0 0 —n m 
The stress-strain relations, then, in x-y system assumes the following 
form, 
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Or Qu Qi Qi 0 0 Cr 
Oy On Q22 Q23 0 0 Ey 
Try } =| Qa Q32 Q3 0 U Yzy (2.46) 
Tyz 0 0 0 Qaa 0 Yyz 
Tzr 0 0 0 0 Qss Yer 
where, 
[9] = [T)* [Q)17] (2.47) 


5. Laminate Analysis 
As discussed earlier, the laminae are stacked to obtain a laminate. In this 
section, the theory associated with the mechanics of laminates is described. 
Consider a laminate composed of N lamina. For the k** lamina of the 


laminate, Equation (2.46) can be written as follows: 


Or Es 
Oy » Ey 

EON (2.48) 
Tyz Yyz 

Tzr K Ver K 


where all matrices must have the subscript K due to orientation of the particular 
lamina with respect to the plate x-y coordinates. 


The functional form of the displacement for a plate are given by: 


u(z,y,z) = u(z,y)+za(z,y) (2.49) 
v(z,y,z) = w(z,y) + zB(z,y) (2.50) 
w(z, y) = Wwo(T, Y) (2.51) 


where u,, v, and w, are middle surface in plane displacements, & and 8 are related 


to the rotations. 


l5 


Substituting Equation (2.49), (2.50) and (2.51) into the Equation (2.24) 


results in: 


Qu, da 




















a "or (2:528) 
E = x y A (2.53) 
ez (2.54) 
fy = > E 4 e) (2.55) 
= ; (3 + w) (2.56) 
Er = > (a + 2 (2.57) 
The mid surface strains are given by the following relations: 
„= = (2.58) 
E, = E (2.59) 
PE > E Hi 2 (2.60) 


The curvature terms, associated with transverse bending are written in 


the form, 


0a 


Kr = Oz (2.61) 
m A (2.62) 
_ 1 / 06 aß 


On substituting the strain-displacement relations into the stress-strain 


relations, we obtain stresses in terms of displacement components, 
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Or Ero + ZKr 


Gy ㅣ A AS 

Try = 이 K) Tro + 2Kgy (2.64) 
tyz Yyz 

Tzr K Yzr K 


For plate/shell type structures, we define stress resultants Nz, N,, N,,, 


Qz, Qy, Mz, M,, and Mzy, as shown in Figure 2.3. For the kin layer, we may write, 


N, Or 
Ny u2 | “y 
ry =] Ory dz (2.65) 
Q. I: 
Q. Ozzr 
For a laminate composed of N-layers, the normal stress resultants are 
given by 
N, N pt, Or 
NE / o, | dz (2.66) 
IN Ey kl tk- Ory 








or, in terms of strains, we have 








NE N tę E Ero tk D Kr 
NS / 이 x | o dz+ 이 x | ~ zdz (2.67) 
3 k=1 | Ytx-1 Ez fk—1 Key 





The equation (2.67) may be written in the form 


[N] = [A] [eo] + [B] [x] (2.68) 


with the elements of the matrices [A] and [B] given by 


N 
A; —ONBOC (to — tg) (5,3 —:1,2,3) (2.69) 
R=l 
Bis = 056 ti) (5 =1,2,3) (2.70) 
k=1 
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The moment resultants are obtained for a K* layer as follows: 


Me t/2 Or 
AS ], A zdz 
My al a y 
or, 


(2:71 
[M] = [B] [eo] + [D] [x] (2.72) 
N 
where D = TS CHR (tk = tka) (2.73) 
D. 


[B], THE STRAIN-DISPLACEMENT MATRIX 


1. Introduction 


In this section, we describe how the matrix [B] may be calculated. 


The matrix [B], which relates the strains and displacements at nodal 
points of an element may be obtained in the form, 


te} = [B] {u} 


(2.74) 
It may be observed that this matrix depends on the choice of the nodal 


degrees of freedom, the shape function and the form of strain-displacement relations. 
In the discussion that follows, we describe the nine noded bi-quadratic 


Lagrangian isoparametric element. There are five degrees of freedom associated with 
each node. 


2. Shape Functions and Their Derivatives 


A finite element is called an isoparametric element if the same interpola- 
tion functions define both the geometry and displacements. 


Figure 2.4 shows the element in the mapped space. For the in-plane 
deformation, the geometry may be interpolated 


18 





Figure 2.3: Positive Directio 
for a Lamina 


ze + 


YT 


ns for Stress Resultants and Stress Couples 
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y 


and the displacements given by 


Ñ |- [N] {z1 y1 T2942 *** T9Y9) 


| 2 = [N] {uy vy uz v2 --- ug vg] 


where the matrix of shape functions is given by 


_[M 0 N, 0 
BE N 0 N, 


We may also write the displacements in a summation notation as 


9 
u = X Nu, 


| 


9 
v = 3. Nivi 


=] 


Ng 0 


0 


Ng 


(2.75) 


(2.76) 


(2.77) 


(2.78) 


(2.79) 


In the case of plate bending, three more degrees of freedom are added to 


the planar displacements. 


The transverse displacement w , 0, and 0, are the rotations of the normal 


to the undeformed middle surface in the z — z and y — z plane, respectively. 


9 
wS 3. Nw 
a 
m = X N, 
t=] 


9 
0, = Y NO, 
1=1 


(2.80) 
(2.81) 


(2.82) 


The shape functions and their derivatives, which are needed for the com- 


putations of strains, are given in Table 2.1. 
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Figure 2.4: Lagrangian Plate Element 


2i 





Figure 2.5: Degrees of Freedom at a Typical Node i. 


It may be noted that the rotational degrees of freedom are treated as 
independent quantities, following Mindlin theory [Refs. 37, 38, 39]. 
3. Jacobian Transformation Matrix 


The chain rules for differentiation of N(&,7) with respect to € and 7 gives, 


ON, 5 ON; Ox ON; Oy 


EN a ^ 
Ot 0729 00 (Se 
ON; ON; Ox ON; Qy 
— = ——— +  — „84 
On Ox On Oy On GR 
or, in matrix form, 

ON]. [Oz dy] pan: 

0 — 1595 OE E" 

OR, | z | oz dz [| 8. esd 

On On On Oy 


The Cartesian derivatives are given by 
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$ a 
| 





Figure 2.6: Displacements in the x-z Plane 
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TABLE 2.1: Shape Function 
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and Their Derivatives 











ON; Ox Oy És ON; 





0t 0 0 
ÊN, |=| ô: a| | aN. 
Oy On On On 


(2.86) 


The Jacobian matrix [J], which contains the derivates of the global coor- 


dinates with respect to the local coordinates is given by, 


oz By 
M< | oz 85 
dn On 


(2.87) 


By using isoparametric element concepts [Ref. 1], we may write the Ja- 


cobian elements as, 





J2 = S 
Tis 
422 = y AD 


The inverse of the Jacobian matrix can be expressed as 


se 1 Jz —J 
J 1 — 22 12 
| | | | | —Ja Jn | 


with the determinant of the Jacobian given by 
| J |= JnuJz2 — Ja Ji? 
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(2.88) 


(2.89) 


(2.90) 


(2.91) 


(2.92) 


(2.93) 


The Cartesian derivatives are, then, given by 


ON, _ 1 
ðr | |J| 


ON; 1 ON, ON, 
s (Ia FE 4 I) 


dy |J] 
4. Formation of [B] Matrix 


ON; 
J 22 DE 


OŚ 


The equation (2.74) may be written as 


Ex 


E, 


2 
ðr 
- | 
o 
Yzy Em 








a 
y 
be 


ON; 
z Ja zę) (2.94) 
5r (2.95) 
| y | (2.96) 


Replacing the displacements in terms of nodal displacements, we get 


al 


or, in matrix form 


i} = [Bi] tu} 


where [B;] is given by 


[Bj2 | 0 
For plate bending, we have, 
Er 0 - 
a i= | uw 
d 
"ry 0 dy 
or, in terms of nodal displacements, 
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8) N; 0 

Ex Or 

z [2lo^ ame 
Ja M 9 
y 


ña | (2.97) 


7 t 
N; 

Or 

(2.98) 


(2.99) 


(2.100) 


eleseo 


OSN; 
Ex 0 — Q 0 | W; 
37 O 0 EE Us (2.101) 
Yzy 0 e NDS Oyi 
y z 
The transverse shears are expressed by 
a 2, N;w; 
| Yyz | = | oy ue | — Y, Nibi (2.102) 
ba e EEEE NA. 
or, 
85 Ni ii L. 
| (> | ES | 8 > N; x aii | 0,; (2.103) 
Z OT m 2 Ni 0 0, 
We may write [B;] for plate elements as 
o <% 0 
Or 
0 0 — 
IB]=|0 -5 —N, (2.104) 
rU —N; 
ae - 0 J 


E. GAUSSIAN QUADRATURE 
1. Introduction 
In this section, we discuss aspects of numerical integration. Here we 
describe the Gaussian method to calculate the integrals, which is widely used in 
finite element work [Refs. 38, 39]. 
2. Summary of Gauss Quadrature 


To approximate an integral J, 


1 
= / ddé (2.105) 
where $(ć) is a function defined in that region. As shown in the Figure 2.7, we 


sample ¢ at the midpoint of the interval and multiply by the length of the interval. 
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Generalization of Equation (2.105) leads to the formula 


T= Wis (2.106) 
I = Wid(€) + Wod(Eo) +--+ Wad (En) (2.107) 


where €; is the location of the integration point z with respect to the origin, W; is 
the weighting factor for point z and n is the number of points at which ¢(€) is to be 
evaluated. 

The sampling points and weights for Gaussian Quadrature, which is 
adapted in the present study, are listed 1n Table 2.2. 


In two dimensions, we may approximate I by 


I - Y. Y WW;d(Es nj (2.108) 
x 
In Eguation (2.6) and (2.7) each coefficient in the integrand matrix must 


be integrated as described above. 
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n 


Figure 2.7: Gaussian Quadrature (€ Coordinate) using Three Sampling 
Points 
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0 774596692 
00 


O 611363116 
O 1109810426 


0 906] 794459 
0.538409 WOT 
00 


U 9374605142 
0 66120932685 
0.2386191861 


0 9491079123 
Q 7415311556 
Q 4088481514 
00 


Q 9602898565 
Q 7966664774 
( 5255324099 
O 1 kK34 140425 
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U RESBRENERY 


QO 247854848] 
Q 6821451539 


0 2369268851 
0 4786086705 
O. SOSARASREO 


0.1713244924 
O 3607615730 
0 46751 39346 


Q 1294849662 
() 2797083915 
0 3818300505 
0 4179591837 


O. 1012285303 
0 22238103415 
O 3] 37066459 
0 3626837834 





TABLE 2.2: Coefficients for Gaussian Quadrature 


III. PROGRAM IMPLEMENTATION 


A. INTRODUCTION 


After having formed the element stiffness matrices, the global matrices are 


formed in a standard manner [Refs. 1, 38], resulting in, for static analysis, 


{F} = [K] {u} 
where {F} : External loads vector 
{K} : Stiffness matrix 
{u} : Nodal displacements vector 
If we know {F} and [K], {u} may be computed for a static problem. Typical 


steps for static analysis may be described as follows: 

Step 1: Input the material properties, plate coordinates, loads, the number of 
integration points, boundary conditions. 

Step 2: For the composite, input number of layers, elements and fiber orien- 
tations. 

Step 3: Determine element matrices in global coordinate system. 

Step 4: Assemble the element matrices. 

Step 5: Solve for displacements and stresses. 


The Figure 3.1 shows the flow chart of the program. 


B. SOLUTION PROCEDURE 

In order to obtain numerical solution, the element matrices were programmed 
and incorporated into an existing finite element program. The implementation in 
double precision was done on a 32-bit Apollo D-3000 series computer. 


Following steps describe the element subroutine. 
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( BEGIN ) 
ㆍ = 






N=1,NO.OF LAYERS 


CALCULATE 


INTEGRATION (2*2) OR (3*3) 


SHAPE FUNCTION 
JACOBIAN MATRIX 


B MATRIX 
FORMULATION OF K 





Figure 3.1: Flow Chart (Main Program) 


Step 1: Material properties from main program are read. 
Step 2: Calculate [Q] according to number of layers. 
Step 3: Select integration point 2 x 2 or 3 x 2 or 3X 3. 
Step 4: Establish shape function. 

Step 5: Determine the Jacobian matrix and [B]. 

Step 6: Formulation of [K]. 

Step 7: For each integration point, do steps 4, 5, and 6, and accumulate [K ]. 
Step 8: Calculate [K] for each element. 
Step 9: Return to main program. 


The flow chart for element level computation is given in Figure 3.2. 
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BEGIN 







N = 1,NO.OF LAYERS 


CALCULATE Q 


INTEGRATION (2*2) OR (3*3) 


* 


SHAPE FUNCTION 
JACOBIAN MATRIX 


B MATRIX 
FORMULATION OF K 


2 


Figure 3.2: Flow Chart (KQML9) 
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IV. NUMERICAL EXAMPLES 


A. INTRODUCTION 

This section discusses the validation of the element formulation and imple- 
mentation by means of selected numerical examples. The isotropic elastic plates are 
solved under different boundary conditions and loadings. 

This is followed by application to selected laminated composite plates to check 
the formulation for such applications. The effects of thickness and integration 


schemes also are investigated. 


B. SIMPLY SUPPORTED SQUARE ISOTROPIC PLATE 
1. Simply supported plate under concentrated load. 
Consider a rectangular isotropic plate under a central point load. The 
geometry and the boundary conditions are shown in the Figure 4.1. The structure 
is modeled using the double symmetry. 


The properties of material / are given by 


E x 10.92 x109 (psi) 
possc 

L = 101m 

P = 400 (lbs) 


The geometric boundary conditions of the simply supported plate are 
w = 0 on all edges. 

As shown in the Figure 4.2, results depicting maximum displacement 
obtained using 2 x 2 integration points, 3 x 3 integration points and ‘heterosis’ [Ref. 
4] elements are compared. As the thickness decreases, the element is more effective 
for integration schemes and for thick plates, the error 1s approximately 2096. It is 


possible that by increasing the number of elements, it may improve the solution. 
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Figure 4.1: A Rectangular Isotropic Plate under Concentrated Point 
Load 


36 


W/W ana! Vs. No. of elements are shown in Figure 4.3. As the number of 
elements are increased, the solution converges to the analytically predicted values. 
It may be noted that reduced order of integration is more flexible, and approaches 
the analytical solution as an upper bound. 

2. Simply Supported Plate under Distributed Load 

Next we consider a rectangular plate under a distributed load. The ge- 
ometry and the boundary conditions are the same as in the previous example. The 
model of the structure, using the symmetry, is shown in Figure 4.4. 

The material properties for this example are given earlier as material I. 

The normalized maximum deflection is shown plotted against the number 
of elements in Figure 4.5. It may be observed that different integration schemes 
converge towards a lower bound. Figure 4.6 depicts the effectiveness of this element 
in predicting behavior of thick and thin plates. 

This element gives better predictions for thick plates than heterosis ele- 


ment, however, for thin plates, heterosis element appears to be better. 


C. CLAMPED-CLAMPED SQUARE ISOTROPIC PLATE 


Consider a rectangular plate clamped on all sides, under a distributed load. 
The geometry and boundary conditions are shown in Figure 4.7. The boundary 
conditions for this problem are implemented by prescribing all degrees of freedom 
to zero, on all four edges. The structure, again, is modeled using the symmetry. 

The material properties for this plate are same as for material I. 

The results showing the normalized maximum displacements are shown in 
Figure 4.8. The heterosis element results in about ten percent error while the 


Lagrangian element gives about twenty to twenty-three percent for thick plates. 
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Figure 4.2: W/Wanat vs. Log (L/T) on Simply Supported Square Plate 
under Central Point Load. 
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Figure 4.4: A Rectangular Isotropic Plate under Distributed Load 
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Figure 4.7: A Clamped-Clamped Rectangular Plate Under Distributed 
Load 


However, for thin plates, the predictions are about the same using either elements. 


Figure 4.9 shows the convergence characteristics of the element. 


D. CANTILEVERED ISOTROPIC PLATE 
Next example considered is a cantilevered plate. The geometry and boundary 
conditions are shown in the Figure 4.10. The material properties, (material II), are 


given by 


E = 1x 14k) (4.1) 
p (4.2) 
L = 10(in) (4.3) 
t = 0.2(in) (4.4) 


The result of this example is shown in the Figure 4.11. As the number of 


elements is increased, the results converge, for all integration schemes. 


E. SIMPLY SUPPORTED LAMINATED PLATE 
1. Graphite-epoxy 
Consider a rectangular composite plate under a sinusoidal load. Two 


types of construction are used. 
Case I : Graphite-Epoxy 0°/90°/90°/0° 
Case II : Graphite-Epoxy 0°/ —60°/60°/0° 


Material properties are given by 


(712 

a < = 4. 
p. - 06 (4.5) 
— = Q. 4. 
z 0.5 (4.6) 
Vi? = 0.25 (4.7) 
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Figure 4.8: 


45 


0» 


oc 


26 ÑZ 


I=ĆIYNY IAA 4 UI 


Če 
ANTOTT 


SLNSIGTT3 AO “ON 
ya 02 91 


zt 





of Element for 


No. 


Nondimensional Deflection vs. 


Clamped-Clamped Plate under Distributed Load 


Figure 4.9: 


46 


0.01(Lbs) 






A 
EDGE 4 | 
^ 
r 0.04(Lbs) 
E y 
DES EDGE 3 
EDGE 2 


Figure 4.10: Cantilevered Plate 
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Figure 4.11: W/Wiana!) VS- No. of Element for Cantilevered Plate 
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či 


gc (4.8) 


The plate is subjected to a transverse sinusoidal load of intensity 


q = qo Sìn E sin T (4.9) 


This problem is studied to compare finite element results with those of 
shear deformation theory [Ref. 24] and elasticity solution [Refs. 16, 17, 18]. The 
results are in Figure 4.12 for Case I. The deflection obtained from the finite element 
solutions agree very well with the exact elasticity solutions. The solutions agree 
well for thin plates, while a large discrepancy between the present predictions and 
Reference 24 is observed. 

The maximum deflection is plotted as the number of elements are in- 
creased in Figure 4.13. The convergence trend may be noted as the number of 
elements are increased. In Figure 4.14, the results for the second type of composite 
(0°/ — 60°/60°/0°) is depicted. 

2. Glass-Epoxy 

This section considers rectangular composite plate under sinusoidal load. 
The two types of construction are Case III, 0°/90°/90°/0 and Case IV, 0? / — 
60°/60°/0°. The plates are square and simply supported. 


The material properties are given by 


Er —— 

p = 06 (4.10) 
(722 

| 8 (4.11) 
M2 = 0.25 (4.12) 
Ey 

m (4.13) 
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Figure 4.12: Central Deflection W,,,, of 4-Layered Square Plate 


(Case I: Graphite-Epoxy) 
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Figure 4.13: Wo, vs. No. of Elements for 4-Layered Square Plate 


(Case I: Graphite-Epoxy) 
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(Case II: Graphite- Epoxy) 
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NO. OF ELEMENTS 


The plate is subjected to a transverse sinusoidal load of intensity. 


= 0, sin sin (4.14) 


Figure 4.15 shows the maximum deflection of this element and of shear 
deformation theory [Ref. 24] and elasticity solution [Refs. 16, 17, 18] for analytical 
solution. Good agreement is obtained for thin plates while discrepancies exist for 


thick plates. Convergence characteristics are shown in Figures 4.16 and 4.17. 
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Figure 4.15: Central Deflection Winar for 4-Layered Square Plate 


(Case III: Glass-Epoxy) 
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Figure 4.16: Wmaz vs. No. of Elements for 4-Layered Square Plate 


(Case III: Glass-Epoxy) 
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Figure 4.17: Wymar vs. No. of Elements for 4-Layered Square Plate 


(Case IV: Glass- Epoxy) 
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V. CONCLUSIONS 


This study is directed towards understanding the linear static analysis of plates 
composed of both isotropic and laminated composites. 

The formulation is based on the principle of virtual work. A finite element 
formulation is presented as a model for the analysis of laminated anisotropic plate 
structures. Several numerical examples for the isotropic elastic plates are solved for 
different boundary conditions and loadings. 

The results show that bi-quadratic isoparametric lagrange plate bending ele- 
ment is effective for analysis of laminated anisotropic plates. Numerical solutions 
agree well with analytical solution. Further, 1t is observed that as the number of 
elements are increased, the convergence to analytical solution is assured. 

The element predicts good results for thin plates while large discrepancies exist 
for thick plates and calls for further investigation. In general, a 3 x 3 integration 
seems to predict the solutions well. 

More numerical experiments need to be done regarding selective integrations, 
and apply to a variety of layer configurations for composite plates. Another area is 
in including the nonlinear terms for buckling problems. The free-vibration analysis 


to predict the natural frequencies and mode shapes is an obvious extension. 
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